iso_boundaries = iso_gdf.boundary
print("Boundary geometry type:", iso_boundaries.geom_type.value_counts())Boundary geometry type: MultiLineString 7
Name: count, dtype: int64
February 11, 2026
matplotlib to layer maps.boundary.centroid.iterrows()zorderBasic map of all ISOs:
We can use spatial methods to get something much nicer-looking. Here’s the end result we’re aiming for:
boundary method returns GeoSeries with the boundary of every geometry in iso_gdfMultiPolygon into MultiLineStringfig, ax = plt.subplots()
iso_boundaries.plot(ax=ax, color="black", linewidth=0.6)
iso_gdf.plot(ax=ax, facecolor = 'blue', alpha = 0.3).set_axis_off()fig, ax so we can plot layers on top of each othermatplotlib aesthetic options like color, alpha = 0.3, linewidth = 0.6 to customize lookDiscussion question: now that we’ve modified the transparency and added boundaries, what observations can you now make about the ISOs?
zorder = 3: places centroids on top of other two layers# previous layers of plot...
# format and place labels
for _, row in iso_gdf.iterrows():
ax.text(x=row.centroid.geometry.x, y=row.centroid.geometry.y, ha="left", va="center",
s=row["NAME"],
fontsize=5, zorder=3,
path_effects=[
pe.Stroke(linewidth=2, foreground="white"), # white halo
pe.Normal(),]
).iterrows(): returns each row’s index + data one-by-one
for _, row: ignores the index, just gets data# previous layers of plot...
# format and place labels
for _, row in iso_gdf.iterrows():
ax.text(x=row.centroid.geometry.x, y=row.centroid.geometry.y, ha="left", va="center",
s=row["NAME"],
fontsize=5, zorder=3,
path_effects=[
pe.Stroke(linewidth=2, foreground="white"), # white halo
pe.Normal(),]
)row.centroid.geometry.x and row.centroid.geometry.y: extract the x and y coords of centroid of that rowha, va: horizontal and vertical alignment relative to x,y# define (fig, ax) as before -- see source code
# format and place labels
for _, row in iso_gdf.iterrows():
ax.text(x=row.centroid.x, y=row.centroid.y, ha="left", va="center",
s=row["NAME"],
fontsize=5, zorder=3,
path_effects=[
pe.Stroke(linewidth=2, foreground="white"), # white halo
pe.Normal(),]
)s: string to use for labelzorder = 3: places labels on top of two other layerspath_effects: adds white stroke on top of labelsax = ax: add to same plot as previous layerszorder = 2: ensures 48 states are plotted below ISOscolor = 'lightgray': color as gray so as to not distract from main focus: ISOsModify the map in spatial2_dps.qmd so that:
firebrickcentroid to place labelsboundary to highlight featuresmatplotlib tools:
zorder to specify order of layersiterrows() to add text labels one by onebuffer.shpgeopandas loads it all in as GeoDataFrame city_fips cbsa_fips city st_fips pop_2010 \
0 0455000 38060 Phoenix 04 1445632.0
1 0477000 46060 Tucson 04 520116.0
2 0644000 31080 Los Angeles 06 3792621.0
3 0666000 41740 San Diego 06 1307402.0
4 0667000 41860 San Francisco 06 805235.0
geometry
0 POINT (4188179.251 -30022.362)
1 POINT (4254222.09 -231751.035)
2 POINT (3618422.756 259366.151)
3 POINT (3676500.155 65351.387)
4 POINT (3454644.566 819013.232)
print("CRS of us_cities_gdf:", us_cities_gdf.crs)
print("CRS of us_states_gdf:", us_states_gdf.crs)
print("CRS of us_powerplant_gdf:", us_powerplant_gdf.crs)CRS of us_cities_gdf: EPSG:3347
CRS of us_states_gdf: EPSG:4326
CRS of us_powerplant_gdf: EPSG:4326
us_cities_gdf does not match us_states_gdf or us_powerplant_gdf!fig, ax = plt.subplots(figsize=(8, 6))
us_states_gdf.boundary.plot(ax=ax, color='lightgray', linewidth=1)
us_cities_gdf.plot(ax=ax, markersize=1.5)
ax.set_axis_off()
plt.show()us_cities_gdf = us_cities_gdf.to_crs(us_states_gdf.crs)
print("CRS of (new) us_powerplant_gdf:", us_powerplant_gdf.crs)
fig, ax = plt.subplots(figsize=(8, 6))
us_states_gdf.boundary.plot(ax=ax, color='lightgray', linewidth=1)
us_cities_gdf.plot(ax=ax, markersize=1.5)
ax.set_axis_off()
plt.show()CRS of (new) us_powerplant_gdf: EPSG:4326
Can we just use sjoin between us_powerplant_gdf and us_cities to find the nearby power plants?
powerplants_near_cities_gdf = gpd.sjoin(us_powerplant_gdf,
us_cities_gdf,
how='inner', predicate='intersects')
print("Shape:", powerplants_near_cities_gdf.geometry.head())GeoSeries([], Name: geometry, dtype: geometry)
No – we get an empty GeoDataFrame!
sjoin hereus_powerplant_gdf and us_cities_gdf are both Point geometries – generally won’t overlap unless they have the exact same coordinatesus_cities_gdf into a Polygon geometry using .buffer()sjoin with us_cities_gdf[Axis(name=Geodetic latitude, abbrev=Lat, direction=north, unit_auth_code=EPSG, unit_code=9122, unit_name=degree), Axis(name=Geodetic longitude, abbrev=Lon, direction=east, unit_auth_code=EPSG, unit_code=9122, unit_name=degree)]
unit_name = degree – not meters!
So we need to re-project us_cities_gdf to a CRS with unit meters. We’ll use EPSG 5070 (NAD 83)
Source: epsg.io/5070
us_cities_100km_gdf = us_cities_gdf.to_crs(epsg=5070)
us_cities_100km_gdf["geometry"] = us_cities_100km_gdf.buffer(100000)
us_cities_100km_gdf = us_cities_100km_gdf.to_crs(epsg=4326)us_cities_gdf into EPSG 5070
us_cities_100km_gdf geometry type is still Pointgeometry with 100km buffer
Polygonus_cities_100km_gdf back into EPSG 4326 to be consistent with other layersus_cities_100km_gdf city_fips cbsa_fips city st_fips pop_2010 \
0 0455000 38060 Phoenix 04 1445632.0
1 0477000 46060 Tucson 04 520116.0
2 0644000 31080 Los Angeles 06 3792621.0
3 0666000 41740 San Diego 06 1307402.0
4 0667000 41860 San Francisco 06 805235.0
geometry
0 POLYGON ((-111.01793 33.71736, -111.00635 33.6...
1 POLYGON ((-109.82395 32.2893, -109.81389 32.20...
2 POLYGON ((-117.34169 34.31756, -117.32295 34.2...
3 POLYGON ((-116.07184 33.02341, -116.0548 32.93...
4 POLYGON ((-121.33763 37.99624, -121.31331 37.9...
us_cities_100km_gdf has inherited the same attributes as us_cities_gdf, except now geometry == POLYGON
bufferfig, ax = plt.subplots(figsize=(8, 6))
us_states_gdf.boundary.plot(ax=ax, color='lightgray', linewidth=0.5)
us_cities_100km_gdf.plot(ax=ax, color='blue', alpha=0.3, zorder = 2)
us_cities_gdf.plot(ax=ax, markersize=1.5, zorder = 3)
ax.set_axis_off()
plt.show()sjoin with 100km bufferpowerplants_within_100km_gdf = gpd.sjoin(us_powerplant_gdf,
us_cities_100km_gdf,
how='inner', predicate='intersects')
print(powerplants_within_100km_gdf.geom_type.value_counts())Point 3325
Name: count, dtype: int64
Discussion questions:
powerplants_within_100km_gdf be based on?sjoinfig, ax = plt.subplots(figsize=(8, 6))
us_states_gdf.boundary.plot(ax=ax, color='lightgray', linewidth=0.5)
us_cities_100km_gdf.plot(ax=ax, color='blue', alpha=0.3, zorder = 2)
powerplants_within_100km_gdf.plot(ax=ax, markersize=1.5, zorder = 3)
ax.set_axis_off()
plt.show()Index(['country', 'country_long', 'name', 'gppd_idnr', 'capacity_mw',
'latitude', 'longitude', 'primary_fuel', 'other_fuel1', 'other_fuel2',
'other_fuel3', 'commissioning_year', 'owner', 'source', 'url',
'geolocation_source', 'wepp_id', 'year_of_capacity_data',
'generation_gwh_2013', 'generation_gwh_2014', 'generation_gwh_2015',
'generation_gwh_2016', 'generation_gwh_2017', 'generation_gwh_2018',
'generation_gwh_2019', 'generation_data_source',
'estimated_generation_gwh_2013', 'estimated_generation_gwh_2014',
'estimated_generation_gwh_2015', 'estimated_generation_gwh_2016',
'estimated_generation_gwh_2017', 'estimated_generation_note_2013',
'estimated_generation_note_2014', 'estimated_generation_note_2015',
'estimated_generation_note_2016', 'estimated_generation_note_2017',
'geometry', 'index_right', 'city_fips', 'cbsa_fips', 'city', 'st_fips',
'pop_2010'],
dtype='object')
After the sjoin, we now have gwh from power plants and city names
generation_by_city = (powerplants_within_100km_gdf
.groupby('city_fips')['generation_gwh_2019']
.sum()
.reset_index(name='total_generation_gwh_100km')
)
us_cities_gdf = us_cities_gdf.merge(
generation_by_city,
on='city_fips',
how='left'
)city_fips, sum up generation_gwh_2019 among nearby power plantspandas merge to merge this new variable back in to us_cities_gdfus_cities_gdf city_fips cbsa_fips city st_fips pop_2010 \
0 0455000 38060 Phoenix 04 1445632.0
1 0477000 46060 Tucson 04 520116.0
2 0644000 31080 Los Angeles 06 3792621.0
3 0666000 41740 San Diego 06 1307402.0
4 0667000 41860 San Francisco 06 805235.0
geometry total_generation_gwh_100km
0 POINT (-112.08906 33.5715) 71285.407290
1 POINT (-110.87812 32.15433) 4824.438325
2 POINT (-118.40647 34.11347) 31489.796582
3 POINT (-117.12247 32.8307) 4202.462840
4 POINT (-122.44305 37.75616) 25909.678010
fig, ax = plt.subplots(figsize=(10, 12))
us_states_gdf.boundary.plot(ax=ax, color='lightgray', linewidth=0.5)
us_cities_gdf.plot(
ax=ax, color='steelblue', edgecolor= 'black', alpha=0.5,
markersize=us_cities_gdf['total_generation_gwh_100km'].fillna(0) / 10,
zorder = 2)
us_cities_gdf.plot(ax=ax, markersize=5, color= 'black', zorder = 3)
ax.set_title("Total Power Generation Within 100km of Large U.S. Cities (GWh, 2019)", fontsize=14)
ax.text(0.5, -0.04,
"Marker size is proportional to total power generation within 100 km (GWh).", transform=ax.transAxes, ha='center', fontsize=10, color='dimgray')
ax.set_axis_off()
plt.show()Encodes size of marker around each city to scale with total_generation_gwh_100km
buffer) to do calculations that don’t show up in the final plotWe’ve covered several use cases for advanced spatial methods:
There are many more methods that we won’t have time to cover.